namespace hms::swe
{

template<typename MeshType>
template<bool singleCell, BlockLocation loc>
void SolverThread<MeshType>::computeBlock( Block<loc>&& block ){
	// std::stringstream s; s
	// << "Thread #" << omp_get_thread_num()
	// << ", time step #" << solver().time().nSteps()
	// << ", block location: " << magic_enum::enum_name(loc)
	// << ", block start: " << block.xBeg() << ", " << block.yBeg()
	// << ", blocksize " << block.nx() << "x" << block.ny()
	// << "\n";
	// std::cout << s.str();
	Timer t;
	const MeshType& mesh { solver().mesh() };
	const ArrayXXs
		& dValues { solver().d().values() },
		& zValues { solver().z().values() },
		& vValues { solver().v().values() };

	/* map local vars to required sizes */
	Array3XsMap q      { block.map(m_qBuf) };
	Array3XsMap fluxes { block.map(m_fluxBuf) };

	auto dBlock { block.template readBlock<1>( dValues ) };
	auto zBlock { block.template readBlock<1>( zValues ) };
	auto vBlock { block.template readBlock<2>( vValues ) };

	ArrayXXsMap
		hBlock { m_hBuf.data(), block.nxRead()  , block.nyRead() },
		qBlock { m_qBuf.data(), block.nxRead()*3, block.nyRead() };

	// #ifndef BLOCKWISE_INTEGRATED_VECTORISE_FLUXES
	// Array1XsMap h { block.map(m_hBuf) };
	// auto d { dBlock.reshaped( 1, block.nRead() ) };
	// auto z { zBlock.reshaped( 1, block.nRead() ) };
	// auto v { vBlock.reshaped( 2, block.nRead() ) };
	// #endif

	/* initialise q, fluxes and h */
	computeStorage(qBlock, hBlock, block, dBlock, zBlock, vBlock);
	fluxes = 0;

	/* initialisation timing */
	setTiming(m_execTimings, LoopStep::init, t);
	// m_execTimings[static_cast<Index>(LoopStep::init)] += t.pop();

	/* compute fluxes */
	#define DECLARE_RESHAPED_VARS \
	Array1XsMap h { block.map(m_hBuf) }; \
	auto d { dBlock.reshaped( 1, block.nRead() ) }; \
	auto z { zBlock.reshaped( 1, block.nRead() ) }; \
	auto v { vBlock.reshaped( 2, block.nRead() ) }; 

	if constexpr ( singleCell ){
		DECLARE_RESHAPED_VARS
		computeFluxes_cell(block, q, fluxes, h, d, z, v, mesh);
	} else {
	#ifdef BLOCKWISE_INTEGRATED_VECTORISE_FLUXES
		computeFluxes_vec(
			mesh, block, q, qBlock, fluxes, dBlock, zBlock, hBlock, vBlock
		);
	#else
		DECLARE_RESHAPED_VARS
		computeFluxes( block, q, fluxes, h, d, z, v, mesh );
	#endif
	}
	#undef DECLARE_RESHAPED_VARS

	/* flux timing */
	setTiming(m_execTimings, LoopStep::fluxes, t);

	/* compute sources */
	addMassSources( block, q );

	/* mass source timing */
	setTiming(m_execTimings, LoopStep::sources, t);

	/* compute friction */
	if ( solver().settings().useFriction )
		computeFriction(
			block, fluxes,
			block.template writeBlock<1>( dValues ),
			block.template writeBlock<2>( vValues )
		);

	/* friction timing */
	setTiming(m_execTimings, LoopStep::friction, t);

	/* calculate new state */
	/* add fluxes (including friction modification) to storage */
	q += fluxes;

	/* get new state from storage */
	auto dBlockNew { block.template writeBlock<1>( this->dNew() ) };
	auto vBlockNew { block.template writeBlock<2>( this->vNew() ) };
	computeNewState(block, q, dBlockNew, vBlockNew);

	/* new state timing */
	setTiming(m_execTimings, LoopStep::newState, t);

	/* calculate maximum allowed time step */
	m_timestep = std::min(
		m_timestep,
		timestep( block, dBlockNew, vBlockNew )
	);

	/* time step timing */
	setTiming(m_execTimings, LoopStep::timestep, t);
}


} // namespace hms::swe